Explore sexual size dimorphism in hummingbird species from Costa Rica
Evaluate Rensch’s Rule
You can have the sections listed here, for instance:
# read data
dat <- as.data.frame(read_excel("./data/raw/supplementary data.xlsx"))
# read tree
tree <- read.tree("./data/raw/imputed hummer tree 339 spp.tre")
dat$`scientific name`[dat$`scientific name` == "Phaetornis guy"] <- "Phaethornis guy"
dat$`scientific name`[dat$`scientific name` == "Phaetornis longirostris"] <- "Phaethornis longirostris"
dat$`scientific name`[dat$`scientific name` == "Heliomaster longisrostris"] <- "Heliomaster longirostris"
dat$scientific_name <- gsub(" ", "_", dat$`scientific name`)
tree$tip.label[tree$tip.label == "Amazilia_amabilis"] <- "Polyerata_amabilis"
tree$tip.label[tree$tip.label == "Hylocharis_eliciae"] <- "Chlorestes_eliciae"
tree$tip.label[tree$tip.label == "Chlorostilbon_canivetii"] <- "Cynanthus_canivetii"
tree$tip.label[tree$tip.label == "Elvira_cupreiceps"] <- "Microchera_cupreiceps"
tree$tip.label[tree$tip.label == "Calliphlox_bryantae"] <- "Philodice_bryantae"
tree$tip.label[tree$tip.label == "Amazilia_edward"] <- "Saucerottia_edward"
tree$tip.label[tree$tip.label == "Amazilia_saucerottei"] <- "Saucerottia_hoffmanni"
tree$tip.label[tree$tip.label == "Amazilia_saucerottei"] <- "Saucerottia_hoffmanni"
tree$tip.label[tree$tip.label == "Amazilia_candida"] <- "Chlorestes_candida"
tree$tip.label[tree$tip.label == "Elvira_chionura"] <- "Microchera_chonura"
# setdiff(unique(dat$scientific_name), tree$tip.label)
subtree <- drop.tip(phy = tree, tip = setdiff(tree$tip.label, unique(dat$scientific_name)))
names(dat) <- gsub(" ", ".", names(dat))
subtree2 <- subtree
subtree2$tip.label <- gsub("_", " ", subtree2$tip.label)
ggtree(subtree2, layout = "circular", col = viridis(10)[7]) + geom_tiplab(size = 4) +
theme(plot.margin = unit(c(30, 10, 30, 10), "mm"))(note that in the data supplied ln.SEX.weight variables are in natural log)
# model settings
v.cv.tree <- ape::vcv.phylo(subtree)
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
prior(student_t(3, 0, 20), "sd"), prior(student_t(3, 0, 20), "sigma"))male_x_mod <- brm(ln.female.weight ~ ln.male.weight + (1 | gr(scientific_name,
cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
prior = prior, iter = iter, control = list(adapt_delta = 0.99,
max_treedepth = 15), file = "./data/processed/male_vs_female_size")| b_prior | sd_prior | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | normal(0, 10) | student_t(3, 0, 20) | 10000 | 4 | 1 | 5000 | 0 | 0 | 10742.36 | 11192.23 | 821034228 |
| Estimate | Rhat | Bulk_ESS | Tail_ESS | CI_low | CI_high | |
|---|---|---|---|---|---|---|
| Intercept | 0.198 | 1 | 10742.36 | 11192.23 | 0.085 | 0.316 |
| ln.male.weight | 0.838 | 1 | 15190.31 | 13900.16 | 0.780 | 0.896 |
female_x_mod <- brm(ln.male.weight ~ ln.female.weight + (1 | gr(scientific_name,
cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
prior = prior, iter = iter, control = list(adapt_delta = 0.99,
max_treedepth = 15), file = "./data/processed/female_vs_male_size.RDS")| b_prior | sd_prior | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | normal(0, 10) | student_t(3, 0, 20) | 10000 | 4 | 1 | 5000 | 0 | 0 | 11919.28 | 13094.22 | 1392122190 |
| Estimate | Rhat | Bulk_ESS | Tail_ESS | CI_low | CI_high | |
|---|---|---|---|---|---|---|
| Intercept | -0.148 | 1 | 11919.28 | 13094.22 | -0.291 | -0.008 |
| ln.female.weight | 1.139 | 1.001 | 15511.00 | 14873.94 | 1.062 | 1.218 |
male_vs_ssd_mod <- brm(aboluteSSD ~ ln.male.weight + (1 | gr(scientific_name,
cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
prior = prior, iter = iter, control = list(adapt_delta = 0.99,
max_treedepth = 15), s = "./data/processed/male_vs_abs_SSD.RDS")| b_prior | sd_prior | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | normal(0, 10) | student_t(3, 0, 20) | 1000 | 4 | 1 | 500 | 0 | 0 | 1486.691 | 1315.158 | 631200404 |
| Estimate | Rhat | Bulk_ESS | Tail_ESS | CI_low | CI_high | |
|---|---|---|---|---|---|---|
| Intercept | 0.009 | 1.001 | 1486.691 | 1315.158 | -0.095 | 0.111 |
| ln.male.weight | 0.065 | 1 | 1871.538 | 1425.567 | 0.010 | 0.121 |
female_vs_ssd_mod <- brm(aboluteSSD ~ ln.female.weight + (1 | gr(scientific_name,
cov = v.cv.tree)), data = dat, family = gaussian(), data2 = list(v.cv.tree = v.cv.tree),
prior = prior, iter = iter, control = list(adapt_delta = 0.99,
max_treedepth = 15), file = "./data/processed/female_vs_abs_SSD.RDS")| b_prior | sd_prior | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | normal(0, 10) | student_t(3, 0, 20) | 10000 | 4 | 1 | 5000 | 1 | 0 | 18472.36 | 14149.97 | 244304415 |
| Estimate | Rhat | Bulk_ESS | Tail_ESS | CI_low | CI_high | |
|---|---|---|---|---|---|---|
| Intercept | 0.053 | 1 | 18472.36 | 14149.97 | -0.068 | 0.168 |
| ln.female.weight | 0.042 | 1.001 | 25317.53 | 14631.72 | -0.027 | 0.112 |
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] viridis_0.6.2 viridisLite_0.4.1 brms_2.18.0
## [4] Rcpp_1.0.9 sketchy_1.0.2 brmsish_1.0.0
## [7] ggtree_3.2.1 readxl_1.3.1 ape_5.6-2
## [10] xaringanExtra_0.7.0 rprojroot_2.0.3 formatR_1.11
## [13] knitr_1.42 kableExtra_1.3.4
##
## loaded via a namespace (and not attached):
## [1] uuid_1.1-0 backports_1.4.1 systemfonts_1.0.4
## [4] plyr_1.8.7 igraph_1.3.5 lazyeval_0.2.2
## [7] splines_4.1.0 crosstalk_1.2.0 ggplot2_3.4.0
## [10] TH.data_1.1-0 rstantools_2.2.0 inline_0.3.19
## [13] digest_0.6.31 yulab.utils_0.0.5 htmltools_0.5.4
## [16] fansi_1.0.3 magrittr_2.0.3 checkmate_2.1.0
## [19] remotes_2.4.2 RcppParallel_5.1.5 matrixStats_0.62.0
## [22] xts_0.12.2 sandwich_3.0-1 svglite_2.1.0
## [25] prettyunits_1.1.1 colorspace_2.0-3 rvest_1.0.3
## [28] ggdist_3.2.0 xfun_0.36 dplyr_1.0.10
## [31] callr_3.7.3 crayon_1.5.2 jsonlite_1.8.4
## [34] lme4_1.1-27.1 survival_3.2-11 zoo_1.8-11
## [37] glue_1.6.2 gtable_0.3.1 emmeans_1.8.1-1
## [40] webshot_0.5.4 distributional_0.3.1 pkgbuild_1.4.0
## [43] rstan_2.21.7 abind_1.4-5 scales_1.2.1
## [46] mvtnorm_1.1-3 DBI_1.1.1 miniUI_0.1.1.1
## [49] xtable_1.8-4 gridGraphics_0.5-1 tidytree_0.4.0
## [52] stats4_4.1.0 StanHeaders_2.21.0-7 DT_0.26
## [55] htmlwidgets_1.5.4 httr_1.4.4 threejs_0.3.3
## [58] posterior_1.3.1 ellipsis_0.3.2 pkgconfig_2.0.3
## [61] loo_2.4.1.9000 farver_2.1.1 sass_0.4.5
## [64] utf8_1.2.2 labeling_0.4.2 ggplotify_0.1.0
## [67] tidyselect_1.2.0 rlang_1.0.6 reshape2_1.4.4
## [70] later_1.3.0 munsell_0.5.0 cellranger_1.1.0
## [73] tools_4.1.0 cachem_1.0.6 cli_3.6.0
## [76] generics_0.1.3 ggridges_0.5.4 evaluate_0.20
## [79] stringr_1.5.0 fastmap_1.1.0 yaml_2.3.7
## [82] processx_3.8.0 purrr_1.0.0 packrat_0.9.0
## [85] pbapply_1.6-0 nlme_3.1-152 projpred_2.0.2
## [88] mime_0.12 aplot_0.1.7 xml2_1.3.3
## [91] compiler_4.1.0 bayesplot_1.9.0 shinythemes_1.2.0
## [94] rstudioapi_0.14 gamm4_0.2-6 treeio_1.18.1
## [97] tibble_3.1.8 bslib_0.4.2 stringi_1.7.12
## [100] highr_0.10 ps_1.7.2 Brobdingnag_1.2-9
## [103] lattice_0.20-44 Matrix_1.5-1 nloptr_1.2.2.2
## [106] markdown_1.3 shinyjs_2.1.0 tensorA_0.36.2
## [109] vctrs_0.5.2 pillar_1.8.1 lifecycle_1.0.3
## [112] jquerylib_0.1.4 bridgesampling_1.1-2 estimability_1.4.1
## [115] cowplot_1.1.1 httpuv_1.6.6 patchwork_1.1.2
## [118] R6_2.5.1 promises_1.2.0.1 gridExtra_2.3
## [121] codetools_0.2-18 boot_1.3-28 colourpicker_1.2.0
## [124] MASS_7.3-54 gtools_3.9.3 assertthat_0.2.1
## [127] withr_2.5.0 shinystan_2.6.0 multcomp_1.4-17
## [130] mgcv_1.8-36 parallel_4.1.0 grid_4.1.0
## [133] ggfun_0.0.7 minqa_1.2.4 tidyr_1.1.3
## [136] coda_0.19-4 rmarkdown_2.20 shiny_1.7.3
## [139] base64enc_0.1-3 dygraphs_1.1.1.6